
#include "GaussLegendre.h"


// Defining possitions in arrays to make things clearer
#define An_0 0
#define An_1 1
#define An_2 2

#define Bn_1 0
#define Bn_2 1

#define Tn_0 0
#define Tn_1 1

#define Pn_0 0
#define Pn_1 1


pthread_t thread[4];

struct GaussLegendre
{
					// CONTENTS:
	pi_type a[3]; 	// a[n-1], 	a[n], 	a[n+1]	| An_0	An_1	An_2
	pi_type b[2]; 	// b[n], 	b[n+1]			| Bn__1	Bn__2
	pi_type t[2]; 	// t[n-1], 	t[n]			| Tn_0	Tn_1
	pi_type p[2]; 	// p[n-1], 	p[n]			| Pn_0	Pn_1
} gl;



void GaussLegendre_A()
{
	// a[n+1] = (a[n] +b[n]) /2
	//gl.a[An_2] = (gl.a[An_1] +gl.b[Bn_1]) /2;
	mpf_add(gl.a[An_2], gl.a[An_1], gl.b[Bn_1]); 
	mpf_div_ui(gl.a[An_2], gl.a[An_2], 2); 
}

void GaussLegendre_B()
{
	// b[n+1] = sqrt(a[n] *b[n])
	//gl.b[Bn_2] = sqrt(gl.a[An_1] *gl.b[Bn_1]);
	mpf_mul(gl.b[Bn_2], gl.a[An_1], gl.b[Bn_1]); 
	mpf_sqrt(gl.b[Bn_2], gl.b[Bn_2]);
}

void GaussLegendre_T()
{
	// t[n+1] = t[n] -p[n] *(a[n] -a[n+1])^2
	//gl.t[Tn_1] = gl.t[Tn_0] -gl.p[Pn_0] *(gl.a[An_0] - gl.a[An_1]) *(gl.a[An_0] - gl.a[An_1]);
	mpf_sub(gl.t[Tn_1], gl.a[An_0], gl.a[An_1]); 
	mpf_mul(gl.t[Tn_1], gl.t[Tn_1], gl.t[Tn_1]); 
	mpf_mul(gl.t[Tn_1], gl.t[Tn_1], gl.p[Pn_0]); 
	mpf_sub(gl.t[Tn_1], gl.t[Tn_0], gl.t[Tn_1]); // t1 = t0-p0*(a0-a1)^2
}

void GaussLegendre_P()
{
	// p[n+1] = p[n] *2
	//gl.p[Pn_1] = gl.p[Pn_0] *2;
	mpf_mul_ui(gl.p[Pn_1], gl.p[Pn_0], 2);
}

void GaussLegendre_Shift()
{
	mpf_set(gl.a[0], gl.a[1]);
	mpf_set(gl.a[1], gl.a[2]);
	mpf_set(gl.b[0], gl.b[1]);
	mpf_set(gl.t[0], gl.t[1]);
	mpf_set(gl.p[0], gl.p[1]);

	/*gl.a[An_0] = gl.a[An_1];
	gl.a[An_1] = gl.a[An_2];
	gl.b[Bn_1] = gl.b[Bn_2];
	gl.t[Tn_0] = gl.t[Tn_1];
	gl.p[Pn_0] = gl.p[Pn_1];*/
}

void GaussLegendre_Init()
{
	// Initializing multiple precision numbers
	mpf_inits(	gl.a[0], gl.a[1], gl.a[2],
				gl.b[0], gl.b[1],
				gl.t[0], gl.t[1],
				gl.p[0], gl.p[1], NULL	);

	// Init [n = 0]
	//gl.a[An_1] = 1;
	mpf_set_ui(gl.a[An_1], 1);
	//gl.b[Bn_1] = sqrt(2)/2;
	mpf_set_ui(gl.b[Bn_1], 2); 
	mpf_sqrt(gl.b[Bn_1], gl.b[Bn_1]); 
	mpf_ui_div(gl.b[Bn_1], 1, gl.b[Bn_1]);

	//gl.t[Tn_1] = 0.25;
	mpf_set_d(gl.t[Tn_1], 0.25);
	//gl.p[Pn_1] = 1;
	mpf_set_ui(gl.p[Pn_1], 1);

	// Calculate first iteration for a, b An_d p
	GaussLegendre_A();
	GaussLegendre_B();

	// Shift all
	GaussLegendre_Shift();
}

void GaussLegendre_Clear()
{
	// Initializing multiple precision numbers
	mpf_clears(	gl.a[0], gl.a[1], gl.a[2],
				gl.b[0], gl.b[1],
				gl.t[0], gl.t[1],
				gl.p[0], gl.p[1], NULL	);
}

void GaussLegendre_CalcPi(pi_type *pi)
{
	//#ifdef TESTING
	//*pi = pow( gl.a[An_1] +gl.b[Bn_1], 2) /(4 *gl.t[Tn_1]);
	//#endif


	// Compute one last t
	GaussLegendre_Shift();
	GaussLegendre_T();

	// Return pi with n+ iterations
	// pi = (a[n] + b[n])^2 / 4t[n]
	//*pi = (gl.a[An_1] +gl.b[Bn_1]) *(gl.a[An_1] +gl.b[Bn_1]) /(4 *gl.t[Tn_1]);
	 mpf_add(*pi, gl.a[An_1], gl.b[Bn_1]); 
	 mpf_mul(*pi, *pi, *pi); 
	 mpf_div(*pi, *pi, gl.t[Tn_1]); 
	 mpf_div_ui(*pi, *pi, 4); // pi = (a1+b1)^2/(4*t1)
}

#ifdef TESTING
void PrintStruct()
{
	gmp_printf("\t \t\t a -> %.10Ff \n", gl.a[An_1]);
	gmp_printf("\t \t\t b -> %.10Ff \n", gl.b[Bn_1]);
	gmp_printf("\t \t\t t -> %.10Ff \n", gl.t[Tn_1]);
	gmp_printf("\t \t\t p -> %.10Ff \n", gl.p[Pn_1]);
	gmp_printf("\n");
}
#endif

void GaussLegendre_Paralelo(long iterations, pi_type *pi)
{
	int i, j;
	void *status[4];

	GaussLegendre_Init();

	for(i=0; i<iterations -1; i++)
	{
		pthread_create(thread +0, NULL, (void * (*)(void *)) GaussLegendre_A, NULL);
		pthread_create(thread +1, NULL, (void * (*)(void *)) GaussLegendre_B, NULL);
		pthread_create(thread +2, NULL, (void * (*)(void *)) GaussLegendre_T, NULL);
		pthread_create(thread +3, NULL, (void * (*)(void *)) GaussLegendre_P, NULL);

		// Syncronizing pipeline
		for(j=0; j<4; j++)
			pthread_join(thread[j], status +j);

		// Debug Output
		#ifdef TESTING
		GaussLegendre_CalcPi(pi);
		gmp_printf("\t %d -> %.10Ff \n", i, *pi);
		PrintStruct();
		getchar();
		#endif
		
		GaussLegendre_Shift();
	}

	GaussLegendre_CalcPi(pi);

	GaussLegendre_Clear();
}

void GaussLegendre_Sequencial(long iterations, pi_type *pi)
{
	int i;
		
	GaussLegendre_Init();

	// Calc pi variables
	for(i=0; i<iterations; i++)
{
		GaussLegendre_A();
		GaussLegendre_B();
		GaussLegendre_T();
		GaussLegendre_P();

		// Debug output
		#ifdef TESTING
		GaussLegendre_CalcPi(pi);
		gmp_printf("\t %d -> %.10Ff \n", i, *pi);
		PrintStruct();
		getchar();
		#endif
		
		GaussLegendre_Shift();
	}

	// Return
	GaussLegendre_CalcPi(pi);

	GaussLegendre_Clear();
}

